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Building on the many existing algorithms for calculating the DC transport properties of quantum 
tight-binding models, we develop a systematic approach that expresses finite frequency observables 
in terms of the stationary Green's function of the system, i.e. the natural output of most DC 
numerical codes. Our framework allows to extend the simulations capabilities of existing codes to a 
large class of observables including, for instance, AC conductance, quantum capacitance, quantum 
pumping, spin pumping or photo-assisted shot noise. The theory is developed within the framework 
of Keldysh formalism and we provide explicit links with the alternative (and equivalent) scattering 
approach. We illustrate the formalism with a study of the AC conductance in a quantum point 
contact and an electronic Mach-Zehnder interferometer in the quantum Hall regime. 
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I. INTRODUCTION 

Numerical simulations of the transport properties of 
quantum tight-biding models are now a mature field 
and are routinely used either alone or in conjunction 
with ab-initio simulations and/or electrostatic simula- 
tions. Examples include a wide variety of systems rang- 
ing from graphene devices^^, quantum Hall effectP^^ 
or spintronic devices^^ to topological insulators^^, 
semi-conducting nanowires^Mj or hybrid superconduct- 
ing systemPHS] 

These simulations are based on two equivalent 
approaches, the scattering (or Landauer-Biittiker) 
formalismP^Sa anc [ the Non Equilibrium Green's Func- 
tion formalisuPH^ (NEGF, itself based on the Keldysh 
formalism). The development of the numerical methods 
for simulating quantum transport is perhaps as old as 
those formalism d 34 " 371 . One of the most popular algo- 
rithms is the Recursive Green's Function methocP^^ 
and its generalizations for tackling larger systems, 
complex geometries and multiterminal measurements 41 . 
Typical observables accessible with those techniques are 
the conductance, the noise, as well as local properties like 
electronic density, current density or spin currents. 

In recent years, it has become clear that probing quan- 
tum systems at finite frequencies can provide entirely 
new physical insights, and finite frequency mesoscopic 
physics is emerging as a field on its own. Among the pi- 
oneer works, the extension of the scattering approach to 
finite frequency initiated irpHEU nas i ec [ to many interest- 
ing predictions including the quantum capacitance, finite 
frequency shot noise, quantum pumping, spin pumping, 
etc. On the experimental level, apart from the huge 
corpus of work with quantum-bits (semi-conductor or 
superconducting based), one observes an increasing in- 
terest in high frequency physics, including several strik- 
ing experi ments on quantum capacitance 4 ^" ^, qua ntum 
inductanc d 52 1 53 1 , or electronic quantum optic^ 4 ^'. 

On the numerical side, finite frequency physics has at- 
tracted comparatively limited interest so far. Pioneer 



works include some calculations by Guo et a pM68J as we p 
as a few otherP^H. 

The aim of this paper is to provide a comprehensive 
approach to numerical simulations of quantum transport 
at finite frequency, bridging the gap between the output 
of the existing numerical codes (typically the Retarded 
Green's function of a sub part of the system) and the 
actual AC observables. We shall provide a systematic 
way to derive formulas that express those AC quantities 
(say the AC conductance) as the trace of a product of 
a few easily accessible Green's functions. Some of these 
formulas are already known in the literature, others are 
new. More specifically, in presence of a perturbation of 
the form V ac cosujt we provide two complementary ap- 
proaches. The first one is a diagrammatic technique for 
a systematic expansion of the retarded Green's function 
in powers of V ac . It allows to easily recover the (mostly 
known) lowest order expressions but also to go beyond 
linear response or the wide band limit using a simple 
set of Feynman like rules. The second method allows to 
access the opposite adiabatic limit (uo — > 0) and system- 
atically expand any observable in powers of uo. 

The paper is organized as follows. Section [TT] contains 
an introduction to numerical simulations of stationary 
(DC) quantum transport. We focus in particular on the 
output of the DC codes that will be the input for the 
calculations of the finite frequency (AC) quantities. We 
are then ready, in Section [ITU to present a rather long list 
of AC expressions that relate one AC quantity (such as 
AC conductance or photocurrent) to the integral (over 
energy) of the trace of a product of DC Green's func- 
tions. For simplicity we give these expressions in the so 
called wide band limit and defer the full expressions (that 
should be used in the numerics) to the Appendix. All the 
expressions given in the Appendix (but the first one) are 
new to the best of our knowledge. In the beginning of 



Section IV we provide some technical details explaining 
how to evaluate the energy integrals in practice. Further 
in this section we consider two applications of increasing 
complexity. In the first one, we consider the AC conduc- 
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tance in a point contact geometry and show that such 
a measurement would allow to extract the time of flight 
of the electrons through the device. In the second, we 
consider the full scale simulation of an electronic Mach- 
Zehnder interferometer made out of the edge states of a 
two dimensional gas in the quantum Hall regime. 

The full machinery used to derive the expressions of 
Section |III| (and the Appendix) is developed in Sections 
|V||VIII| It provides a set of simple rules allowing to de- 



rive the expressions given in Section III as well as any 
other similar expression beyond those given explicitly 
here. Section [V] briefly introduces the necessary nota- 
tions and main results of the Keldysh formalism. Section 
[VI] focuses on the particular case where the perturbation 
is periodic in time. We then proceed with developing 
systematic perturbative expansions around two distinct 
limits: first we provide a diagrammatic technique to ex- 
pand the results in powers of the amplitude of the per- 
turbation. Second, we expand around the adiabatic limit 
where the frequency of the AC perturbation is ve ry small. 
These expansions are presented in Section VII ^pertur- 
bation applied inside the system) and Sectior jVIII (per- 
turbation in the electrodes). 



II. 



BUILDING BLOCK FOR AC TRANSPORT: 
THE GREEN'S FUNCTIONS OF THE 
STATIONARY PROBLEM 




FIG. 1: Sketch of a multi-terminal mesocopic system con- 
nected to M = 4 leads. Numbers with a bar on top denote 
the parts of our system, i.e. for the device region and 1...4 
for the leads. 

We begin this manuscript with a brief account of the 
numerical methods used for simulating the stationary 
(DC) transport properties of quantum systems. We focus 
on tight-binding models within the NEGF formalism^ 
which has become, perhaps, the standard approach for 
this problem thanks to the de velopme nt of recursive^®' 
and other advanced algorithm d 41 * 79 * 8Q 4 The basic objects 
obtained as the output of those numerical tools are the re- 
tarded Green's functions of the system which will be the 
input of the expressions for the AC observables derived 
later in this manuscript. 



We consider a generic quadratic discrete Hamiltonian 
of an open system written as 



(1) 



where (c n ) are the usual creation (destruction) opera- 
tor on site n. The site index n is very generic and includes 
all the degrees of freedom present in the system (spatial, 
momentum, spin, electron/hole, orbital). In terms of the 
(infinite) H matrix, the retarded Green's function is de- 
fined as, 



g(E) = (E + ie-ny 1 , 



(2) 



where e is an infinitely small positive number. We now 
focus on the particular geometry of our mesoscopic sys- 
tem: a finite region, referred as connected to M semi- 
infinite electrodes 1...M, see Fig. [I] In the rest of this 
manuscript, notations such as Afhrn' always refer to the 
corresponding sub-block of the full infinite matrix A. 
Each electrode is kept at equilibrium with chemical po- 
tential /i m and temperature T m . The first step in the 
NEGF formalism consists in integrating out the elec- 
trodes degrees of freedom. One obtains for the retarded 
Green function inside the device region Qqq(E): 



E-H 



M \ _1 

m=l / 



(3) 



where H = Hqo is the Hamiltonian matrix projected in- 
side the device region (see Fig. [I]) and Y. r (m;E) is the 
(retarded) self-energy due to the presence of the lead m. 
The latter is given by, 



Z r (m;E)=H- ^(E + ie-H^)-'H fh - . (4) 

Equations ([3| and Q a re the ty pical raw output of, 
say, recursive techniques^^^ 9 -^ from which one can 
calculate various physical observables such as conduc- 
tance or current noise. For instance, the celebrated Lan- 
dauer formula for t he current flowing from lead m reads 
in this context 29 ^, 



r M 

= \ \ dE (f™ ~ fm') Tr [0or m /£jr. 



(5) 



where we have introduced the Fermi function f m (E) = 
1/ '^(E-^/kTm _|_ the standard broadening matrix 



T m (E)=i[T r (m;E)-T r \m;E) 
and the shortcut 

g l (E) = g m (E+—). 



(6) 



(7) 



In what follows, our aim is to extend expressions such 
as Eq.([5| to AC quantities. Those expressions will involve 
photon absorption and emission processes and therefore 
will require the calculation of Qi{E) for I ^ 0. 
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III. COOKBOOK FOR CALCULATING AC 
OBSERVABLES 

This section is intended as some sort of dictionary 
where we provide expressions for various finite frequency 
response functions. Those expressions allow to extend 
directly DC numerical codes to the AC regime. They 
are g iven without derivation and we refer to Sections \V\ - 



VIII| for the proofs and/or the rules for deriving other 
observables not included in this section. 

The generic perturbation we consider takes the form 



W = coscjt^ W 



(8) 



where W is a rather arbitrary (hermitian) matrix to be 
specified below. This form of perturbation is suitable for 
AC electric field, AC magnetic field or lattice deforma- 
tion. 



A. External AC perturbation 

We first consider the situation where a uniform AC 
voltage is applied to one of the contacting electrodes, say 
contact m! . The perturbation takes the form 



W 



(9) 



where 1^/ is the identity matrix inside region fa' . The 
expressions below are given in the so-called wide-band 
limit (WBL) where the energy dependance of the elec- 
trodes (retarded self-energies) is neglected. The most 
general (but also cumbersome) expressions without this 
approximation are given in Appendix [Aj 

The first important observable is the current I m (t) 
flowing through contact m. It can be decomposed ac- 
cording to the different harmonics of the perturbation, 



I m (t) = TteJ2l m (luj)e- 1 



loot 



(10) 



1=0 



The AC conductance matrix 81 (in the absence of DC 
bias) readPl 



-iS n 



dl m (lu) _ 
dV ac 



dETi 



f(E + M 



(ii) 



ftjUJ 



At small frequency, it simply reduces to the Landauer 
formula 

In the WBL, we find an absence of rectification: 



dVl 



0. 



(12) 



However this is not the case in general, as seen in Eq.( A4) 



in Appendix [AJ Let us briefly discuss what the above re- 
sult implies on the modeling of the electrodes, i.e. on the 



effect of the position of the dashed line in FigjT] (where 
the electrodes end and the system starts). Suppose that 
a real system is made of a small mesoscopic region con- 
nected to contacts. The contacts themselves are made 
of two parts: a very wide metallic part (energy inde- 
pendent, connected to the measuring apparatus) followed 
by a narrower region (energy dependent) which is itself 
connected to the mesoscopic region. From the modeling 
point of view, one must decide where the electrodes start 
(in the wide or narrow region). At the quantum mechani- 
cal level, the position of the electrodes is totally arbitrary 
and simply corresponds to which degrees of freedom are 
integrated out. Therefore at this level, the physics is un- 
affected by the position where one chooses to define the 
electrodes (in the wid e region or in the narrower one). 
The fact that Eq.(A4) gives a non zero result (lead in 
the narrower regi on) o r a vanishing one [lead in the wide 
region where Eq.(A4) reduces to Eq.(12)] therefore in- 
dicates that the difference between the two cases takes 
place at the statistical physics level, i.e. upon assuming 
that each electrode always remains at its thermal equi- 
librium. The correct choice between the two above men- 
tioned possibilities depends on the inelastic mean free 
path: when almost no inelastic collisions take place in 
the narrower region, the electrodes should be considered 
to be in the wide region and the rectification effect van- 
ishes. At higher temperature, the inelastic mean free 
path decreases, and the narrower region eventually be- 
comes thermalized, which leads to a non zero rectification 
effect. 

Another interesting limit is the adiabatic limit when 
the frequency uj is very small while the amplitude of the 
perturbation V ac can remain arbitrary large. To zeroth 
order in Hlj, the current (for m ^ m!) is simply given by 
a trivial extension of the DC result, 



dE 



df_ 
dE 



Tr 



Go^m'Gl^ 



eV ac cos UJt. 
(13) 



This expression is l inea r in V ac in the WBL. However, 
in general (see Eq.( |A5| )) adiabatic current contains all 
higher orders in amplitude as well. More interestingly, 
the first correction to adiabaticity reads, 



teuj 

47T 



j Mi- 



di 
dE 



Tr 



9$, 
' dE 



' dE 



(14) 



eVnr smut. 



Note that while the adiabatic current follows exactly the 
slow changes of voltage, this correction is out of phase. 

Another important observable is the electronic density 
n(i,t) = (c\(t)ci(t)) on site i whose decomposition in 
harmonics reads 



n(i, t) 



(i) + Re ^ lw, m)e 



ilujt 



(15) 



1=0 
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where n eq (i) is the stationnary density in the absence of 
the time dependent potential. We refer to the response 
function dn{i, lco>, m')/dV ac as the generalized injectivity. 
It is a straightforward generalization of the injectivity 
defined ir l 43 | 4T | at small frequency and it reads 76 , 



dn(i, lco>, m') 

dVa C 



£ J dE (f(E) - f(E + M) [g 2 r m ,gl 

(16) 



B. Internal AC perturbation 

Let us now turn to the case where the AC perturbation 
is applied inside the device region, 



W 



(17) 



where the matrix block Wqo can take an arbitrary form 
inside the device region, allowing one to include many 
types of perturbation such as electric gates or time de- 
pendent magnetic field. 

In the WBL the linear in V ac current response is^ 



dl m (luj) ie 2 
dVac h 



dE(f(E)-f(E + hw))Tr V m Q 2 WQl 

(18)' 



The DC (rectified) current contains only even orders 
in V ac (because sign of V ac is just a phase shift in Eq.(|8|). 
Thus, the leading order contribution (in the absence of 
DC bias) readsPES 



d 2 I m {0u 



(f(E) - f(E + M) Tr \g wg 2 TglwglT 



- (f(E -hu)- f(E)) Tr \g wg- 2 Tgl 2 wglT. 



(19) 



where 1 = J2 m T m- 

Similarly, the 2 nd harmonics of the current is given by, 



! dE { 



d 2 I m (2uj) _ ie 3 
dV a % ~ 2h 

(f(E)-f(E + fuv))Tv 

+ (f(E -two)- /(E)) Tr g 2 wg wgl 2 r. 



g 2 wglwgi 2 r r , 



(20) 



A particularly interesting case is the current gener- 
ated upon perturbating the onsite potential on site i (i.e. 
Wfcj = eVaSikSii). The A C response function is the (gen- 
eralized) emissivit ^ 43 * 47 * and it reads^, 



dl m (luj) ie 2 
dVu h 



(21) 



Note that the emissivity defined by Biittiker has a mean- 
ing of the number of particles (rather than the current) 
being emitted into the lead m as a consequence of an 
internal perturbation. Thus, to come back to the origi- 
nal definition one has to multiply Eq.(21) on both sides 
by l/(ieu). Then its relation to the introduced above 



injectivity, Eq.(16), will become obvious. 



Last, we introduce the frequency-dependent Lindhard 
functiorP^H that relates a change of density on site j to 
the perturbation on site j' as^, 



n(iw,j,/) 



dn(j, Iuj) 

dVyy 



ie 
'2tt 



J dE{(f(E)-f(E + hw))[g 2 ] jr 

f(E + m [gl] ; [gl] ; - f(E) \g 2 \. y [g ].,.}, 



Gl] , (22) 

x 3 3 



where we implied an expansion similar to Eq.(|T5|) for the 
electronic density on site j. We note that this expression 
does not assume the WBL. 



C. Coulomb interactions: screening 

The expressions given so far do not take electron- 
electron interaction into account. Besides missing po- 
tentially relevant correlated physics, the non interacting 
expressions violate two important laws: conservation of 
current (some finite AC charge may pile up in the device 
region) and "gauge invariance" (raising the AC voltages 
of all leads simultaneously may produce some finite AC 
current), 



(23) 



However, at a scale large enough, some screening will 
eventually take place to restore the global electric neu- 
trality of the system (and these two laws). A minimum 
treatment of electron-electron interactions therefore im- 
plies solving the transport equations together with the 
Poisson equation, following the general framework devel- 
oped in the group of Buttiker^ 3 -^. In practice, one has 
to solve the Poisson equation for the AC electric poten- 
tial V(r,t) (we use spatial notations r instead of general 
coordinates i in this subsection) which for the first AC 
harmonics reads, 



dV(r, luo^m') 

' dVac = 

I- J dr f U(loj,r : r f 



dn(r, lco>, m r ) 



dV ac 
dV{r',luj,m') 



(24) 



(e: dielectric constant). In a second step, one calculates 
the full AC current response in the leads as, 



dlmO-u) 

dVac. 



— T 

— -L r 



/ 



dr 



dV(r,lw,m') 



dV{r) 



dV a , 



(25) 
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We therefore find that the full knowledge of the injec- 
tivity E q.fL6| ), emissivity Eq.(21) and Lindhard func- 
tion Eq.([22[) i s ne eded to complete a full self-consistent 
calculati^tffilZZI 



D. Relation to the scattering matrix formalism 

The Green's function approach taken here can be 
equivalently recast into the scattering approach devel- 
opped h J 44 l 46 l 81 J. In particular, the scattering matrix 
S(E) is connected to the retarded Green's function Q{E) 
through the Fisher-Lee relatiorP^^. In the WBL, the 
equivalent expressions in the scattering formalism are 
simply obtained by using the formal substitution, 



r n (?ov r m -iS nm (E) - il 



(26) 



where S nm (E) is the scattering submatrix between the 
electrodes m and n. For instance, for the AC conduc- 
tance, Eq.dTTl), the substitution provides, 



p 2 r jrp 

r m , m ,(w) = -- J — (f(E) - f(E + Hlj)) 
x Tr[<5 mt7 



S J mm ,(E)S mm ,(E + hw)], (27) 



which is precisely the expression derived in Ref. 81 . Note 
that the trace in Eq.(27) is taken over the open channels 
in the lead. 

In general (not in the WBL) mapping between the 
two formalisms requires solving the scattering problem in 
presence of the oscillating field via the Floquet approach. 
For more details on this approach and its relation to the 
NEGF see Ref P. 




FIG. 2: (Color online): Real (left panels) and imaginary (right 
panels) parts of the integrands appearing in Eq.( 28 ) as a func- 
tion of energy close to the band edge E = 0. The frequency 
was chosen to be ftw = 0.1. There are pronounced peaks at 
energies E = and E = —fvuj. Upper panels: A§i. Lower 
panels: A21 (black solid line) and A%i (red dashed line). 



with 



^21 



Tr \A^ r (E; E + fvuj)G 2 K r 1 a {E + twj] E)Gi 
A21 = Tr [A r 2 r (E; E + huj)G 2 A r 1 r (E + Hu; E)G ] 
Tr \A% a (E; E + fiw)^A? (£7 + fiw; E)G ] Q 



a aa 
^21 



where 



(29) 
(30) 
(31) 



IV. APPLICATIONS 

In this section, we apply the formalism on three prac- 
tical examples: the AC conductance of a simple one di- 
mensional chain, a quantum point contact (QPC) and 
an electronic Mach-Zehnder (MZ) interferometer in the 
quantum Hall regime. These examples, of increasing 
complexity, are chosen to illustrate how the numerical 
calculations can be performed in practice and how the AC 
physics can provide insights absent in DC. The station- 
ary Green's functions Gi at the root of the AC expressions 
were obtained with the knitting algorithm described in 
RefP. 

The AC observable we consider is the AC conductance, 



which is given by Eq.( Al ). Eq.( Al ) can be written as the 
sum of three terms, 



/ 



dE 



h J fajj 
-f(E)A r 2 r 1 +f(E + hw)Ar i } 



Tr[(f(E)-f(E + hu,))A 



(28) 



At(E;E') = T c (m;E)-T b (m;E'). 



(32) 



(here c, b G {a, r} stands for the retarded or advanced self 
energy). The three terms (29), (30) and (31) are direct 
outputs of recursive Green's function like techniques so 
that the main numerical difficulty lies in the evaluation 
of the integral over energies. 



A. Technical details on the numerical integration 

We start with the AC conductance of a simple one 
dimensional chain described by the Hamiltonian, 



H = 2 



00 



h.c. 



(33) 



(the constant 2 serves to offset the bottom of the band to 
E = 0). The device region is of size L so that we suppose 
that the system stays in thermal equilibrium for n < 
(left lead, region 1) and n > L (right lead, region 2). 
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For a single site L = 1 device, the (onsite) Green's 
function of the system can be easily obtained 85 , 



S(E) 



, forE < 0, 
= ,for|£| <4, 



i v /4-(£7-2)2 

, 1 , for£ > 4. 



(34) 



We find that it contains square root singularities 
1/y/E-Eo which appear on the edges of the band (or 
more generally in quasi ID system, whenever there is 
an opening/closing of a new conducting channel). Typ- 
ical plots of the integrands are shown in Fig. [2] These 
singularities are integrable but may require a very fine 
discretization mesh. In practice, we find that advanced 
integration routines, such as QUAD which is used in this 
work 86 (based on an adaptative mesh and Gauss quadra- 
ture formula), can handle these singularities properly and 
provide precise results within a reasonable computing 
time. A faster convergence is often obtained when one 
provides the routine with a precise location of the po- 
sition of the singularities. This location can be found 
by simple calculations of the DC transmission of perfect 
wires as a function of energy using a dichotomy algorithm 
(the singularities occur at the energy where the transmis- 
sion has a step like increase). An alternative route is to 
remove the singularities using a local change of variable. 
One starts by locating the singularities and dividing the 
integration range in small segments with (at most) one 
singularity on one boundary of the segment. The integra- 
tion dEf(E) on one segment (all those integrations 
can be done in parallel) is then performed using the lo- 
cal change of variable E = y/E — Eq which removes the 

singularity, i.e. one integrates J Q a E ° 2dEg(E) with the 
non diverging function g(E) = Ef(E 2 + Eo) (See also^ 
and the straightforward improvement to include two sin- 
gularities in one segment). 

The AC conductance of the one dimensional wire is 
shown in Fig. [3] for two lengths L = 40 and L = 200. 
We find that the calculation performed keeping only the 
"WBL-like" term (this is the only surviving term 
in the WBL) in Eq.(28), is equivalent to using Eq.(27), 
derived within the scattering approach 8 -^, in the large L 
limit. In order to compare the two approaches, we inte- 
grated numerically (assuming zero temperature) Eq.(|27|) 
using the actual dispersion relation of the Hamiltonian 
( [33] ), E(k) = 2 — 2cos/c, and S21 = exp(z/cL) (per- 
fect transmission). To understand why both formalisms 
coincide when L — >> 00, note that terms A%[ and A™ 
typically oscillate as exp [zbz (k(E + two) + k(E)) L] oc 
exp [±2ikpL\ (kp: Fermi momentum) so their integrals 
quickly vanish when kpL ^> 1 and only the (WBL) term 
A21 oc exp [i (k(E + huj) — k(E)) L] remains (in agree- 
ment with Refs! 81 * 88 ! where the fast oscillating terms are 
neglected in the derivation of the current operator). 




0.05 0.1 
flG) 



0.15 



FIG. 3: (Color online): Absolute value of the AC left-to- 
right conductance for the one-dimensional wire of length 
L — 40 (black cirles) and L — 200 (red rectangles). Sym- 
bols: numerical calculation with the Green's function formal- 
ism, Eq.(28) (kee ping only the A%[ term), lines: scattering 
approach Eq.(27) (see text for details). We find a visible dif- 
ference between the two approaches for the small size L = 40 
that disappears when L increases. Fermi energy is Ef — 0.17. 



B. Quantum Point Contact 

We continue with a quasi one-dimensional wire of 
width W and length L connected to two reservoirs. AC 
bias is applied to the source lead (S) and we are inter- 
ested in the current response in the drain (D), see Fig. 4] 
The Hamiltonian is the direct extension of Eq. ( 33 ) to the 
quasi ID geometry. 



00 w 



H — 2 ^ c I+l,m C n,m + c I,m+l c n,ra + h.C. 



n= — oo m=l 



(35) 



The dispersion relation for this discrete model is (in units 
of the hopping constant) 

E n (k) = e n + 2 - 2 cos k, (36) 

wich corresponds, in the continuum limit, k —> 0, to 



E n (k) 



(37) 



where the transverse energy e n = 2 — 
2 cos (nir/(W + 1)) , n = 1 . . . W, and k is the lon- 
gitudinal momentum. The integer n defines the 
quantized values of the transverse momentum and 
thereby enumerates the conducting channels. 

We focus on the regime where only the first (n = 1) 
channel is open and use the corresponding transverse en- 
ergy ei as our reference energy. In this limit the system 
is described by a unique velocity, hence the time of flight 
through the system is well defined. In the scattering 
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FIG. 4: (Color online): Rescaled amplitude and phase (up- 
per inset) of the two terminal AC conductance Tds 



IT 



Ds\e 



i<j>AC 



for a quasi one-dimensional wire of width W 



10. Lower inset: schematic picture of the setup. Fermi energy 
is chosen to have only one propagating mode. The parame- 
ters are: E F = 0.68e 2 , L = 100 (black circles), E F = 0.79e 2 , 
L = 150 (blue rectangles), Ef — 0.79e 2 , L = 200 (green tri- 
angles), where e 2 is the energy of the second mode opening 
(see Eq.(36)). Different symbols correspond to the numerical 



integration of Eq.( 28), while the lines are calculated with the 
help of Eq.(38) exploiting the full tight-binding dispersion re- 
lation (36 ). Red dashed line is the analytical fit using Eq.(39 ). 



All the lengths are in units of the lattice constant. 



matrix approach, the system is described by its trans- 
mission matrix Sds(E) = ex.p(ikL): a wave packet is 
entirely transmitted and only acquires a (energy depen- 
dent) phase. Then, equation (27) reads 



Yds 



h JE F -huj 



dE 



■ exp {iL [k(E + Huj) - k(E)}} . 

(38) 

In the continuum limit, when hu <^ Ep it can be further 



simplified using Eq.(37) into 



-DS 



( hu 2 r \ 
V 4E F J 



(39) 



where r = L/vp is the time of flight from the source lead 
to the drain (see lower inset on Fig. [4]). From this simple 
calculation we notice that the AC conductance gives ac- 
cess to two characteristic time (or energy) scales of the 
system. Indeed, the numerical results of Fig. [4] indicate 
that the absolute value of the AC conductance and its 
phase can be very well described by this simple scaling 
law up to moderately large frequencies. The scaling pa- 
rameters arising in this case allow for the extraction of 
the time of flight and the longitudinal Fermi energy. We 
note that, as the Fermi energy Ep and thus the velocity 
vp that enter the previous expression are counted from 
ei (i.e. they are in fact the longitudinal Fermi energy 



and velocity), one can actually slow down the electrons 
to bring these times and energy scales into an experimen- 
tally accessible window (GHz range). 

A practical way to implement an effective quasi-one 
dimensional wire is through a quantum point contact 
(QPC) formed by confining an electron gas with elec- 
trostatic gates placed on top of a semiconducting het- 
erostucture. We add an electric potential to our quasi 
one dimensional wire Hamiltonian, 



oo W 



V = V g ^®x( n - n 0)®y( m - m 0) C i,m C n 



n— — oo m=l 



(40) 



where we have used the following confining "saddle 
point" potential (see the lower inset of Fig. [5|i for a color 
plot of the potential), 



<S> y (m) = 



tanh 



n + Vx 



2— tanh 



m + riy 



tanh 



n-rjx 



-tanh 



m-rjy 



(41) 



(42) 



where the parameters £ x ,€y control the steepness of the 
potential (they are chosen big enough so that the po- 
tential rises essentially in an adiabatic way), and (#o,2/o) 
determines the position of the center of the QPC. The 
effective length (width) of the QPC is 2rj x (2rj y ). In this 
case the dispersion relation in the gated region may be 
considered similar to Eq.([38| except that now the trans- 
verse energy e n is determined not only by the width, but 
also depends on the parameters of the QPC (TJ y ^ y ,V g , 
etc.). In our calculations we controlled the value of V g 
(keeping other parameters fixed) to have only one open 
channel (at a given Ep) in the gated region (see upper 
inset on Fig. [5^l) . The results of the numerical simula- 
tions of the AC conductance for this system are given in 
Fig. [5j Fig. [5^i shows the absolute value of the AC con- 
ductance as a function of the driving frequency. Different 
symbols correspond to the numerical results [Eq.(j28|] for 
different values of the gate voltage V g , hence for different 
effective longitudinal velocities, as shown in the upper 
inset. The fitting lines are obtained from the scattering 
matrix formalism, Eq.(38), making use of the dispersion 
relation (|36|). Again, the transverse energy of the open 



mode, €1(77^, V g , ...), was chosen as the energy refer- 
ence. Comparing the two approaches, we see that the 
closer we are to the edge of the propagating mode (vari- 
ous symbols on the upper inset of Fig. [5^i), the worse is 
the fit given by Eq.(38). This is due to the fact that the 



scattering matrix formula is applicable when the Fermi 
velocity is a smooth slowly varying function of energjl^, 
which breaks down near the band edge (Ep w ei). 

From the slope of the phase of the AC conductance 
((f) = ujt, curves similar to the inset of Fig.|4j not shown) 
we can extract the effective time of flight r of the elec- 
trons through the QPC, see Figs. |5b,c. We find that r 
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with a schematic of the two interfering edge states. The 
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FIG. 5: (Color online): Quantum Point Contact, (a) Absolute 
value of the AC conductance as a function of driving frequency 
for different values of V g depicted on the upper inset. Upper 
inset: dc conductance as a function of V g . The lower inset 
represents actual potential profile in the device forming the 
QPC. (b),(c) Extracted time of flight as a function of QPC 
length, 2r/ x , (see Eq.(41)) and V g , respectively. Fermi energy 
in the calculations was Ef = 1.2 providing 4 open channels in 
the leads, while only one being transmitted through the QPC. 
Other parameters: W = 11, L = 124, £ x — rj x /2, ^ y — 3rj y /2, 
and rjy — 2. 



scales linearly with the QPC length 2r] x (ballistic trans- 
port) while it increases as we tune V g towards the closing 
of the propagating mode (the velocity vanishes when the 
mode becomes evanescent, E F < ei). The above calcu- 
lations did not take screening into account. However, as 
all the calculations were performed on perfectly transmit- 
ting channels there is no "Landauer dipole" in the prob- 
lem. Calculations done assuming perfect screening (en- 
forcing current and gauge conservatio n 144 1 89 1 suggest that 
the above physical signatures remain observable in pres- 
ence of screening, and in particular should allow for the 
measurement of the time of flight. 



C. Mach-Zehnder interferometer 

We close this section with a discussion of the AC re- 
sponse of an electronic MZ interferometer in the quan- 
tum Hall regime^HSl The setup consists of the two- 
dimensional electron gas confined in a finite geometry, 
connected to three reservoirs: source (S), drain (D) and 
internal drain (D'). Fig. [6^1 presents the sample together 




0.02 0.04 

fico 



0.01 0.02 

fico 



FIG. 6: (Color online): Mach-Zender interferometer, (a) Car- 
toon of the system with the interfering paths represented 
by the solid and dashed white lines. Length of the sample 
L = 80, width of the system W = W h + 22, where W h is the 
width of the hole. The system is in the quantum Hall regime 
at filling factor v — 1. The two QPCs were chosen to be 
semi-transparent Ti,2 = \- (b) Transmission characteristics 
of the QPCs. The blue dot corresponds to the Fermi energy 
at which QPC is half transparent, (c) AB oscillations of AC 
conductance as a function of magnetic flux through the hole. 
Black circles, red rectangles, and blue diamonds correspond 
to the driving frequency huo = 0.003, 0.007, 0.01 respectively. 
Width of the hole Wh = 45. (d) wd as a func tion of frequency 
for different values of Wh extracted with Eq.(50) (see text for 
details). From bottom up Wh = 30,35,40,45. (e) Function 
|S(a;)|, see Eq.(52), as a function of frequency. Black circles 



and red diamonds correspond to Wh — 35,45 respectively, 
(f) Respective phase of the AC conductance <j)AC — ojt as 
a function of frequency. In plots (c) and (e) all the symbols 
were calculated with the Green's function formalism, Eq.(28), 
modeling the QPCs with Eq.(41). The connecting lines are 



the corresponding analytical fits, Eq.(50), with parameters 



calculated from the Green's function-based numerics. 
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corresponding Hamiltonian of the system is 

n, m£S 

+V n (43) 

where is a Peierls phase^ in the Landau gauge (mag- 
netic flux through a lattice plaquette in units of e/H) 
and S the list of sites in the system (orange region in 
Fig. [6|i plus the semi infinite leads). The system is put 
in a strong magnetic field driving the system into the 
quantum Hall regime. We focused on the case with a 
filling factor v = 1 (single edge channel). The additional 
electrode (D') is necessary to avoid the paths to make 
multiple loops (that would turn the interferometer into 
a Fabry- Perot). The electric potential V n ^ m defines the 
two QPCs in the lower arm of the interferometer (see the 
previous section). The two QPCs play the role of beam 
splitters for the quantum Hall edge state. The inter- 
fering paths are shown by the solid and dashed lines on 
Fig. |6^l. We put an additional Aharonov-Bohm (AB) flux 
<Pab through the hole of the interferometer that allows 
to change the relative phase between the paths without 
changing the edge states. We calculated the AC conduc- 
tance as a function of the AB flux (fAB and frequency u 
as shown in Fig. [6] A similar setup was studied in DC in 
RefP. 

We start with simple analytical considerations using 
the scattering matrix approach for the edge states. Let 
us assume, for simplicity, that both QPCs are charac- 
terized by energy independent transmission probabili- 
ties Ti ? 2 (and corresponding reflection probabilities are 
R12 = 1 — ^1,2) • We describe the QPCs by their scatter- 
ing matrices, which can be parametrized as follows 



Rk 

% 



1% 

/Rk 



1,2. 



(44) 



The source-to-drain transmission amplitude Sds then 
consists of the contributions from two paths, path a (solid 
line on Fig. |6^i) which is a consequence of two sequential 
transmissions through the QPCs and path b (dashed line) 
arising from sequential reflections 



(45) 



Traversing either path, an electron acquires a phase ip a ,b- 
The phase itself contains two contributions, the AB phase 
caused by the magnetic flux through the hole and a dy- 
namical phase of the propagating plain wave along the 
path, 

^ X (E) = k(E)L x + y x ,x = {a, 6}, (46) 

^PAB =Pb-(Pa, (47) 

where k is the longitudinal wave number of the edge state 
and L x is a length of the corresponding path. In order 
to calculate the AC conductance we need to specify the 
energy dependence of the phase factors in Eq.(45). We 



note that by varying energy we modify only the dynam- 
ical part of the phase, while the AB flux is unaffected. 
Thus, we make the following ansatz 



dk 



A 



L X ^(E F ) 



{a,b}. 



(48) 
(49) 



Actual value of i/j x (Ef) and X x depends on the boundary 
conditions defining the geometry of the setup. Assuming 
we know the edge state dispersion relation, we can relate 
the latter to the group velocity of the edge state v g at 
the Fermi level via X x = L x /(Hv g ). 

At this stage we are prepared to calculate the AC con- 
ductance with equations (|27|), ([45]), and (48). The scat- 



tering formula is valid when huj <C Ep and we will use it 
to carry out the integration in energy. Let's assume, for 
simplicity, that the QPCs are tuned to half transparency 
^1,2 = 1/2. Then, after a straightforward calculation we 
obtain 



Y ds (u,<Pab) = ^e Zi 



cos 



sin 



2 



urd 

2 cos(MEf)-MEf)) 
Lb — L a 



Lb + L a 
T = — , , # 



, (50) 
(51) 



We have two time scales naturally appeared, namely the 
average time of flight r and the relative time coming 
from the asymmetry between the paths. 

Now we turn to our numerical results^see Fig. [6| and 
compare them to the simple formula (50). On Fig. [6b 



we plot the DC characteristics of the QPCs considered 
in our modelling. The Fermi level is fixed at a half 
transparency value. We remind that in order to obtain 
Eq.(50) we assumed that the transmission characteris- 
tics of the QPCs are energy independent. However, as 
one can see from Fig. [6}d, in our sample detuning from 
the assumed value becomes important for hu c± 0.01 
and is of the order A7\ 5 2 — 0.15. Next, Fig. [6h rep- 
resents the plots of the AB oscillations of the AC con- 
ductance as a function of flux for three different val- 
ues of driving frequency. Symbols of different types on 
the Figure represent the Green's function based calcu- 
lation, Eq.([28|), while the solid lines are corresponding 
fits according to Eq.([50|. This fit is obtained as follows: 
(i) we perform a DC calculation (which corresponds to 
uj — » in Eq.(50)) and find </>o = arccos(T^5'(0, tt) — 
Yds^O, 0)), which corresponds to the phase offset at 
zero flux cpAB = 0, see Fig.Qc; (ii) we compute the 
phase uo'd = 2 arccos \Yps(w, 7r) + Yds{w, 0)1 a t small lo; 



(hi) finally, we plot Eq.(50) using the extracted d and 
i/jb^Ep) — ^ a (£^i?) = cj)o + (fAB- We see that this formula 
describes quite well the numerical data, especially at low 
frequencies. However when the driving frequency is in- 
creased deviation of the numerical data from the fit be- 
comes more pronounced. A more detailed analysis shows 
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that these deviations can be accounted for by including 
the quadratic term in energy, which we have neglected in 
Eq.([48|). In Fig. [Gjd we plot the extracted ujil) as a func- 
tion of driving frequency. Various plots correspond to 
the samples with a different width of the hole in the in- 
terferometer. Varying this width, we modify the length 
of the upper path (dashed line in Fig. . Extracting 
the corresponding slopes, i.e. for each sample, we are 
able to calculate the velocity of the edge state v g owing 
to Eq.(51 ). For the parameters chosen in our calculation, 
B ~ 20T and the lattice constant ao — lnm, we have ob- 
tained the velocity v g ~ 0.7ao£ft _1 (which corresponds to 
v g ~ 10 6 m/s), where t = h 2 /(2ma%) is a hopping param- 
eter of the tight-binding model 29 used to simulate the 
setup, m- effective mass. 

In Figs. [6^,f we present the AC conductance calcula- 
tions as a function of frequency. On Fig. [6^, using the 
source-to-drain conductance ^(w,^), we plot the 
function 



T ds (uj,tt) - T ds (uj,0) 



r DS (o,TT)-r DS (o,o) 



(52) 



which for a simple case of Eq.([50j) reduces to 
sin(u;$/2)/(a;$/2). Again, the symbols correspond to the 
Green's function calculations, while the lines are given 
by the analytical fit, Eq.([50]), using the calculated be- 
fore values of uo'd (see Fig. |6]i). Black circles and red 
diamonds represent calculations with various values of 
the hole width (allowing to change the length difference 
between the paths). We notice that for the frequencies 
huj < 0.02, the numerical data is very well fit by Eq.([50]), 
while at higher values of hu this is no longer true. There 
are two reasons for this. First, at high enough frequencies 
detuning of the QPCs from the half transparency value 
becomes significant (see Fig. J6j)) and we cannot neglect 
the energy dependence of the transmission/reflection am- 
plitudes. Second, due to the dispersion of the edge state, 
there is always a contribution from the quadratic term 
neglected in the expansion (48), which becomes impor- 
tant at high frequencies. 

Finally, we plot the frequency dependence of the phase 
of the AC conductance (see Fig. |6]f) varying the hole 
width in order to extract the second time scale, r ac- 



cording to Eq.(50). We find again the value of the group 
velocity v g ~ 0.72ao£ft -1 , which is consistent with the 
previous result. 

To conclude, we have analyzed the AC response of the 
MZ electronic interferometer and found that the non- 
equilibrium dynamics makes it possible to address the 
internal time scales of the setup via transport measure- 
ments. 



[A] as well as systematic tools to derive other observ- 
ables. The formalism developed in the manuscript is 
based on the Non Equilibrium Green's Function formal- 
ism (NEGF). While the application of NEGF to sta- 
tionary quantum transport is now standard, we quickly 
review below th e main features of its time-dependent 
version^- 5 * 68 * 70 * 71 !. Our starting point is a general time- 
dependent quadratic Hamiltonian for our system: 



H(t) = ^H nm (t) C t Cm . 



(53) 



We do not include electron-electron interactions besides 
some mean field treatment as was discussed in Sec. IIII CI 
The basic objects that will be manipulated in this pa- 
per are various sorts of Green's Functions (GFs). The re- 
tarded 0; m (*,f), advanced 0£ m (t,f), lesser ®nm(M') 
and greater (5> m (t,t f ) Green's functions are defined as, 



® r nm(t,t') = -j_6(t - t')({c n (t),cl(t')}) , (54) 
Km(t,t') = \e{t' - t)({c n (t),cl(t')}) , (55) 



<5< m (t,t>) = l (cl(t>)c n (t)}, 
®> m (t,t') = - l -(c n (t)cl(t')). 



(56) 
(57) 



where c^it) (c n (t)) corresponds to (c n ) in the Heisen- 
berg representation. 

The retarded and lesser /greater GFs satisfy the follow- 
ing equations of motiorP*E21 



(58) 



ih^ - H(t) ) <S K (t, t') = 0, k=<,> . (59) 



Table [J summarizes the various Green functions intro- 
duced so far, as well as the one needed later for the AC 
formalism. 



A. Dyson equation 



It is often convenient to split the full Hamiltonian ( 53 ) 



into a sum of an "unperturbed part" H and a perturba- 
tion V(t): 



H = ft + V(t). 



(60) 



r. TIME DEPENDENT NEGF FORMALISM IN 
A NUT SHELL 

The rest of this manuscript is devoted to the deriva- 



tion of the expressions given in Section III and Appendix 



This splitting can (and will) be done in several different 
ways, dictated by convenience. For instance, V can be a 
time-dependent potential, or a hopping element between 
the device and the leads, or a sum of the previous two, 
etc. Introducing # r (t,t') and ^ < (t,t / ), the unperturbed 
Green's functions associated to one can derive the 
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TABLE I: Summary of notations. 



Type of GF a 



Description 



K (£, t') GF of the system when the leads and the 



scattering region are decoupled, see Eq.(60) 



g K (t, t') GF of the device sub-block for a system with 
decoupled leads. g K (t,t') = 0g Q (t, t') 

<5 K (t,t') 



Hamiltonian ( 


53 


)• 


G K (t,t?) 


GF of the device sub-block, see Eq.( 


68). 



G K (t,t') = &l- (t,t') 



Gi(E) GF of the device sub-bloc k wi th / photons 
emitted/absorbed, see Eq.(81). 



G^ (E) n-th order in V ac of the device sub-block GF 



with / photons emitted/absorbed, see Eq.(90) 



Qi (E) Retarded equilibrium GF of the system at 



energy E + 



orium tor£ or tne syste 
, see Eqs.ji]) and 



a K = r, a, <, > 



Dyson equationl ^ 33 ' 78 1 , which relate the full GFs to the 
unperturbed ones. They read, 

® r =0 r + r *V*0 r , (61) 

r =0 r + r *V*0 r , (62) 

® K = Q K + Q r * V * <5 K + % K * V * a , k =<, > , (63) 

<5 K = $ K + <5 r * V * % K + K * V * a , k =<, > , (64) 

where the symbol * stands for convolution with respect to 
time and matrix product with respect to the site indices: 

(A * B)ij{t,t r ) = J2f dt A ik {t,t)B kj {t\t r ) (65) 

k 

and V should be understood as S(t — t')V(t) in a convo- 
lution. 



B. Integrating out the electrodes 

From now on we restrict V(t) to the matrix elements 
that couple the leads to the system plus (possibly) a time 
dependent potential in the device region. Then from 
Eqs.(61) and (63) one arrives at 



G r = g r + g r * (S r + V) * G r 



(66) 



and 



G K = Q K + g r * V * G K + g r * S r * G K + g r * E K * G a 
+ <f * V * G a + g K * S a * G a , k =<, > , (67) 

where we have introduced special notations for the (00) 
device sub-block (see Section |TT|) , 

H = H m , V = Voo, G K = eg , 5 K = g o , ... (68) 



and also the self-energies defined as 

M 

E K = ^ S K (m), « = r, a, <, > . 



(69) 



Utilizing Eq.(58), equation (66) can be rewritten in terms 
of an effective equation of motion, 

(ih^ - H^j G r (t, t') - (£ r + V) * G r = S(t - 0, (70) 



while the lesser and greater GFs (67) with the help of 
Eqs.(58)-(59) and (70) are simplified to, 



G K = G r * E K * G a , «=<,>. 



(71) 



In the absence of time-dependent perturbations in the 
Hamiltonian (60), the Fourrier transform of the self- 



energies Yi < (m) and E > (m) are given by fluctuation- 
dissipation theorem which results ir J 29 * 32 * 33 ^, 

T < (m;E)=if m (E)T rn (E), (72) 
Z>(m]E) = -i(l-f m (E))T m (E), (73) 

with T m defined in Eq.Q. 

C. Expression for the current 

The current associated with the ra-th lead is found us- 
ing the approach Q f 31 l 7Q l 71 lZ8] ; i. e> calculating the change 
of the number of particles in the lead due to connection 
with the device region and thereby with other leads too. 
So, the expression for the current reads 

I m = -e<^>, N m =J2 cUt)c*(t). (74) 



Taking into account the definition (56), Eq.(74) trans- 
forms into 



I m (t) = eJ2 (V ai (t)<8< a (t,t) - V ia (t)®ai(t,t)) . (75) 

i<Gcr , 

Let us introduce auxiliary quantity, 



J(m;t,t') =eJ2 {V ai (t')<8<(t, t') - V iQ (i)C5< (t, t')) 

(76) 



Exploiting Eqs.(|63j) and (64), one can get 

6^0 = tin™ * V^o * ©00 + Q^fn * VfnO * <Sg 



(77) 



®L = Via * V 0m * Qrirn + «50 * V 0m * Cm • (78) 
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Substituting expressions (77) and (78) into Eq.(76), we 
come to 



J(m; t, t') = eTr (G r * S< (m) + G< * S a (m) 
-E r (m)*G < -E < (m)*G a ). 



(79) 



Having this, we can easily find our expression for the 
current, since 



Im(t) J(m\t,i). 



(80) 



VI. EFFECT OF A PERIODIC POTENTIAL 

We now use explicitly the fact that the perturbation ([sj) 
is periodic in time. We introduce the Wigner coordinates 
{r = t - t',T = (t + t')/2} and notice that the GF is 
periodic in T, G r (r,T) G r (r,T + 2tt/uj). Thus it 
is possible to expand the GF into a Fourier series with 
respect to T and into a Fourier integral with respect to 



J68I73I9QI 

Gr(T ' T) = ISh £ e~i ET e~ iujlT Gl(E). (81) 

Z= — oo 

We will use extensivley the fact that when C(t, t') — A*B 
is a convolution of two quantities, one gets^, 

Cl {E)= E AiAE+^)B l2 (E-^). (82) 

h+h=l 

For instance, the currents reads, 
dE 



Im(t) 



/^E ^JKm;^, (83) 



l— — oo 



where using Eq.(82) for each term in Eq.(79) we arrive 
at 



Ji(m;E) = e ^ Tr 



G r h (E-\ 



-^)G<(E-^)-E<(m;E- 



(84) 



This is the main starting point of all the subsequent 
derivations. We still have to supplement it by three equa- 
tions. First, using Eqs.(71) and (82) it is straightforward 
to find 



M 



m' = l l\+l2+lz=l 

Ki^E-^ + ^GUE-^-^). (85) 



VII. PERTURBATION IN THE DEVICE 
REGION 

In this section we consider the case when the pertur- 
bation is applied inside the scattering region. Therefore 
we assume that the leads are at local thermal equilibrium 
and unaffected by the perturbation. This implies that V 
in Eq.(70) is given by Eqs.([8|) and (17). Performing the 



transformation (81) in the Eq.(70) we get 



Second, if we consider the definition of the advanced and 
the retarded GFs and a symmetry of the transformation 
(81 ), we easily get 



G?(E) = [G r _ l (E)] i . 



(86) 



Finally, with the two previous expressions and Eq.(84) 
one can deduce 



Ji(m;E) = [J_«(m;£7)]T. 



(87) 



In the next sections we develop a systematic way to cal- 
culate the current (83), i.e. to express the GF elements 



G\(E) in terms of known quantities. We consider two 
type of perturbations: internal perturbation (in the de- 
vice region) and external perturbation (in the contacts). 



E 



hujl 
~Y 
w 



H — Y. r (E 



hujl 



eV ac -GU(E ■ 



or equivalently, 



~2' 



cm 

w 



eV ac --G\+ x {E 



~2 > 



G\{E)-eV ac g l {E)\(G\_ l {E- h ^) 



+GWi(E- 



fibJ . 



SloQo(E). 



(89) 



In the two following paragraphs we explore two comple- 
mentary limits: small perturbation amplitude (eV ac ) and 
adiabatic limit (uo — >• 0). 
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A. Limit of a small perturbation amplitude 

Let us assume that the perturbation amplitude is much 
smaller than all characteristic energy scales in the system, 
e.g. hopping constant between sites. Then we can solve 
Eq.(88) iteratively in powers of eV ac . The solution takes 
form 



(90) 



n=0 

Gl n \E)=g l (E)^(Gl n _- 1 \E 



(91) 



for n > 1 and 



Gf\E) = 5 lfi g l {E). 



(92) 



This equation can be solved iteratively. It is instructive 
to write down explicitly the first and second order con- 
tributions, 



G^{E) 



ME^Q-xiE) + S^G-iiE^GiiE), 



(93) 



W 



w 



G?\E) =6 l , 2 G2(E)^Go(E)^G-2(E) 

w w 

+5i,-2G-2{E)—G (E)—G2{E) 

( w w 

+S lfi lGo(E) Y G- 2 (E) Y Go(E) 



w w 

G {E)-G2(E)-Go{E) 



(94) 



The structure of the two previous equations suggests the 
following simple diagrammatic representation of an ar- 
bitrary order contribution (see Fig. [7]). The diagrams 
are made by horizontal propagating lines [Qi{E)\ sepa- 
rated by "photon absorption/mission" vertical wavy lines 

[W/2]. In order to build G z (n) (E) one has to remember 
the following "Feynman rules", 

• To get the contributions of order n, draw n wavy 
lines pointing up or down in all possible configura- 
tions (there are 2 n diagrams). 

• Each diagram of order n gives a contribution to 

(E) where I is the difference between the num- 
ber of up wavy lines and down wavy lines. 

• Read the diagram from left to right. Starting from 
Qi{E), each wavy line correspond to a factor W/2 
followed by another Qv{E) with V decreased by 2 
(up wavy line, a "photon" is emitted) or increased 
by 2 (down wavy line, "photon" absorbed). Repeat 
until the end of the diagram. 



1-st order 

+i 



+i 

r 



2-nd order: 

+2 



2 0+2 



to , of 



-2 



+2 



FIG. 7: First and second order diagrams. Numbers over the 
straight lines correspond to the index / for Qi. Direction of 
the wavy line tells us whether the absorbtion or emission of 
a photon takes place. Each diagram corresponds directly to 
one term in Eqs.( 93|94 ) 



B. Adiabatic limit and beyond 

We now turn to a very different limit where the driving 
frequency (and not the amplitude) is the small parame- 
ter of the problem. When the perturbation varies slowly 
the system follows adiabatically. We introduce the gen- 
erating function F as, 



F(z,E)= zl °^ E 



l= — co 



Hu)l N 



(95) 



A closed equation for F can be obtained by expanding 
the self energy, 

TT{E + huX) = Z r (E) + ^r^T- ( 96 ) 

Assuming that we work around the wide band limit and 
using Eq.(88) (evaluated at energy E + foujl/2), we obtain 
up to first order in the derivative of Y. r (E), 



F + hujzF ad 



dF 



= F n 



with 



E - H -Vr(E) - eV ac f (z + \)' 



(97) 



(98) 



Note that F a d corresponds to the adiabatic limit: when 
one evaluates F for a given time T (see Eq.(81)), 
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Fad(e~ lUjT i E) corresponds to the stationary retarded 
Green's function at energy E for the potential at time T, 
i.e. assuming that at a given time T, the potential varies 
so slowly that it can be considered as constant. Higher 
order terms can be obtained straightforwardly and corre- 
spond to higher derivatives of F. For instance, to second 
order, one should add the following to the left hand side 
ofEq.d97) 



(M 2 F , ^ 
2 ad[ ] dE 2 



> d 2 F 
dz 2 



OF 

dz 



(99) 



Equation (97) allows for a systematic calculation of F 
(and therefore the G/), for instance by expanding it in 
powers of Huj. To first order we get, 



drn dFad 
dE\ dT ' 
(100) 

And higher orders are obtained straightforwardly. We 
emphasize that in the adiabatic limit the processes con- 
tain a (arbitrary) large number of absorbed/emitted pho- 
tons, hence the role of the generating functions which 
correspond to an instantaneous basis. The resulting ob- 



servables are given in the cookbook section III and in 
Appendix [Al 



where % is the Hamiltonian of the leads and device (when 
they are decoupled) and, 



M 



sin u)t f 



c]c a ) (103) 



is the new coupling between them. By doing this change 
of basis we are back to the situation where the leads are 
kept at (local) thermal equilibrium, whereas the effect 
of the perturbation is completely transferred to the cou- 
pling matrix between the latter and the scattering region. 
As a result, the self energies of lead fa' now acquire an 
additional phase factor, 



E"(ra';M') = ^{m r -t-t f )e- % -^ 



1 (sin ojt — sin cut') 



, K = r,a,<, 
(104) 



where E^ra'; t — t') is the equilibrium self-energy (in the 
absence of the AC field). Expanding the phase factor in 
Eq.(104) in terms of Bessel functions, 



+oo 



n= — oo 



VIII. PERTURBATION IN THE LEADS 

The formalism developed above can be extended to 
homogeneous perturbations in the leads. The algebra is 
very similar with one notable exception: multiple absorp- 
tion/emission processes are now allowed. We suppose 
(for definitness) that a bias voltage V ac cos out is applied 
to lead fh' , see Eq.(|9|. 

A. Equation of motion 

It is convenient to change the basis in the lead affected 
by the perturbation in order to move to a frame where 
the lead is stationary. The AC voltage then gives rise 
to a time-dependent phase factor in the coupling matrix 
between the lead and the device. This is easily seen with 
the help of the unitary transformation 



U = exp 



eV ac cos(ut')N 



,fi= E c « c - 



(101) 

The Hamiltonian after the transformation refers to the 
old one, 



the tra nsform ation (81) applied to the (perturbed) self- 
energy (104) gives 



£f(m';25) 



71= — OO 



(105) 



Finally, the equation of motion has the form 



M 



HujL 



E ^(m'-E+^il-h^GU^E 



Huj 1 1 



l 1= —oo 



(106) 



This equation is the starting point for the approximation 
schemes considered below. 



B. Limit of a small perturbation amplitude 



H' = UHU* -ihU 



Then, it consists of 



H'(t)=H + V(t), 



(102) 



Let us now expand G\(E) in powers of eV ac /huj <C 1, 
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The first o rder r esult can be obtained by a direct expan- 
sion of Eq.(fT05l), 



Ef(m' ] E) = Z K (m' ] E)5 lt0 



+ 



eVac 
2?vjj 

eV ac 



^(m'-E-^-T^m'-E- 



~2 



2hw 

+ 0(e 2 V* c ), K = r,a,<, 



(108) 



where we used the power series representation of Bessel 
functions. Utilizing the notation introduced in Eq.(32), 
we obtain 



A3 A" As 

+4 < < +4 < -4 
Tl T2 tl 

FIG. 8: Example of a diagram that contributes to G^ 2 \E). 
Each diagram is characterized by the set of upper and lower 
numbers which correspond respectively to the number of emit- 
ted and absorbed photons. The numbers along the horizontal 
line are calculated from the diagrammatic rules. 



G<°>(25) 
G?\E) 



^,000, 



(109) 



)Q- 



(110) 



which is similar to expressions ( 92 ) and ( 93 ) . 



In analogy with the case of internal perturbations, a 
systematic diagrammatic expansion can be constructed. 
The set of rules to obtain all the contributions to (E) 
is given by, 

• Draw all diagrams with 1 < p < n wavy lines. 
Contrary to the previous case, the wavy lines now 
point both up and down, reflecting the possibility of 
multiple absorption and emission processes. Each 
wavy line i is associated with two positive integers 
and n\ (nf-\-nf > 1) that correspond to the two 
types of processes. We have 



i=l 



(in) 

(112) 



2 = 1 



• Read the diagram from left to right. Starting from 
Qi (E) , each wavy line corresponds to a factor T l n a n& 
(see below) followed by another Qv(E) with V = 
/ + 2{n ( * — nf ) (nf "photons" are emitted and nf 
are absorbed). Repeat until the end of the diagram. 



• The "vertex" is given by 



1 



1 



xAl^ IE 



-[l 



(113) 



with the matrix A^, (E) given by 
Al,(£) = £(-!)* 

i=0 

xA£, (E + (2i - q)^; E + (2i + 2-q) 2 




(114) 



(A 



where I y I is the binomial coefficient. As an example, 
Fig. ^corresponds to one contribution to the 12-th order 

( 12) 

with I = 4, namely G\ (E). Using the above rules, this 
diagram gives, 

g[ 12) (e) = g 4 Tl 3 g T%g 4 Tl 5 g. 4 , (115) 



where, according to Eq.(113), 
1 



rr-i4 

^1,3 — 



2 4 3! 



[AZ,(E-fiu;E)-3AZ,(E;E + ^) 



+3A£, (E + hou;E + 2hu) - A£, ( E + E + 3fiw)] , 

(116) 

T °n = xtL [A£, (E\ E + huj)- KZ> (E + hw;E + 2huj)} , 

(117) 



^2.0 - 2 2 2 ! 

1 



- 5A^, (E - 2Hu- E-hw) + 10A^ (E - Hoo; E) 

- 10A£, (E; E + hJ) + 5A£, (E + hw;E + 2hw) 

(E + 2hcu; E + 3hu)} . (118) 

We see that in contrast to the case when the pertur- 
bation was inside the scattering region, multiple-photon 
absorption/emission processes are allowed. This fact can 
be also understood from the concept of the sidebands 
(with energy shifted with respect to the Fermi level by 
an amount ±nftw) which have been i ntroduced in the 
context of AC scattering theor}^ 9 -* 73 * 91 1 . 
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C. Adiabatic limit and beyond 



Finally, we consider the adiabatic limit following a sim- 



ilar procedure to the one presented in section VII B The 



procedure is very similar except for the expansion of the 
self energy of the lea d und er AC perturbation. The gen- 
eral expansion of Eq.(105) reads, 



m':E 



(l-h) 



^ k\ 

k=0 



d k Y. r (m , ;E) 
dE* 



{l-2l 1 -2n) k J h+n 



eVac 



J n 



eVac 



(119) 



r 



which, using the following two identities for Bessel func- 
tions, 



Jn+l(x)J n {x) = S it0 , (120) 
J n (x) = J n -l(x) + J n -l(x), (121) 



n= — oo 

2n 



allows to obtain the self energy to any value of k. Re- 
stricting to first order, we get 

E^m'; E + ^(Z - h)j « i Ilj0 r(m'; E) 

eV ™fx ■ x ,dT r (m f ;E) fkul <9Z r (ra'; E) 
TTWiA + °h-i) + °h,o — 



dE 



This expansion corresponds to the wide-band limit (k 



Q; b6 | 7U | and its first correction (k — 1). It is expected 
to be very accurate in metallic leads, for instance. At 
this level, we introduce again the generating function for 
F(z), Eq.(95) and obtain the same equation 97 as for the 



internal perturbation case provided one replace W by 
[— dY. r (m f \ E)/dE]. The corresponding results can hence 



be adapted to this case straightforwardly. Note that be- 
yond this first order closed equation can also be obtained, 
but this simple replacement rule does not apply anymore. 



IX. CONCLUSIONS 

Numerical simulations of quantum transport has be- 
come an ubiquitous tool for mesoscopic physics and are 
more and more commonly used to help the design of na- 
noelectronic devices. On the other hand there is a general 
trend of mesoscopic physics and microelectronics towards 
GHz or even higher frequencies, so that developing a gen- 
eral framework to tackle finite frequency transport is be- 
coming of increasing importance. In this manuscript we 
have developed the corresponding formalism allowing to 
derive a large set of formulas that express AC observ- 
ables in terms of numerically accessible quantities. We 
provide systematic rules to construct other expressions 
that we did not give explicitly. Hence, this manuscript 
can either be used as a recipe book for extending DC 
numerical tools to AC, or as a starting point for further 
developments. 
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Appendix A: General results for AC observables 
given in Section |III| 

In this Appendix we present the formulas for AC ob- 



we will split them in two parts: response to the internal 
and external perturbations. 



servables beyond the WBL, see Section III Just as before 



1. External perturbation 



First, we consider the situation when the AC voltage is 
applied to one of the contacts, say m! . Below we present 
various system response functions to such a perturbation 
assuming there is no DC bias in the system. 

The AC conductance matrix is given by 



T m yM = ^ J ^l^f(E)-f(E + hLo))TT[A%(E;E + ^)g 2 A™,(E + hw;E)gi 
-f(E)Tv [AZ(E; E + fiu)g 2 A!Z,(E + huj; E)g ] + f(E + huj)Tv [\%(E; E + tko)g\AZ{E + huj; E)Q\ 

— Tr [f(E)D 1 - f(E + huj)D 2 ] <5 m>m , 1 



(Al) 



Functions Di^ in the diagonal part are denned as 

D l =i(p 2 - 9l) r(m'; E) + (q - gfj AZ(E; E + hu), 

(A2) 



The stationary (DC) component of the current to lead- 
ing order in perturbation amplitude reads (m ^ m') 



D 2 =i(g 2 - ad) r(m'; E + hco) 
+ (g 2 -gl)A%(E;E + tko). 



(A3) 



d 2 I m (0u 
dV a % 



le 



dE 



-g AZ,(E; E + tiw)g 2 A™,(E + friu; E)GlT m (E) 
(f(E -fiu)- f(E)) Tr [g AZ,(E; E - foj)0_ 2 A£, (E - fko; E)QlT m {E) 
-g A%,{E;E- hw)gl 2 AZ(E - hw;E)glT m (E) 
-if(E)Tv [g (T m ,(E + hw) + T m ,(E -tiu)- 2T m ,{E))glT m {E)\ } 



(A4) 



Adiabatic current. Assuming that the frequency of the perturbation is very small, to zeroth order in ftjjj the cur- 



19 



rent is (m ^ m') 



f(t) = lj dE{f(E)Tr [i ^°\t,E)-F^\t,E))T m (E)-F^(t,E)T(E)F^(t,E)T m (E)] 



+Tr 



Fi0)(tjE) ( ^M f{E) + %T m ,(E)) F^(t,E)T. 



V dE 



dE 



eV nr COS Ut 



}■ 



(A5) 



where the definition of the functions F^'^\t,E) are Correction to the adiabatic current. It is th e leading order 
given by, correction of the order ~ fiw to Eq.(A5). Note that it 

F^(t,E) 



contains all orders in V n , 



E-H - E r (E) ■ 



9Tr< ^' E) eV ac coswi' 
(A6) 



F^(t,E) = iF^ 



dE \ dE 



dE / dE 



2 



(A7) 



<(*) = SlKt) + 8I 2 m (t) + 5I 3 m (t) + 5I 4 m (t) + 5I 5 m (t), 

(A8) 



where separate parts read, 



5I 1 m {t) = e J j/(£)Tr [i (f™ (t, E) - F^ (t, E)) T m (E) - F« (t, E)T(E)F™ (t, E)T m {E) 
F^\t,E)T{E)F^\t,E)T m {E) 



F^(t,E) (^E^lf {E ) + ^r ra> (£0) F^(t,E)T m (E) (A9) 
+F(°\t,E) ^M m + ^T m ,(E)j F^(t,E)T m (E) 



si 2 m (t) 



eh f dE , 
2ttH' 



I 



-Tr 



d(F«»(t,E) + F«»*(t,E)) ( dT m {E) df 
dt \ ~~dE~^ ' + dE Tm{E) 



(A10) 



si: 



-Tr 



eh f 

ar m >(E) 



eh f dE d 
2Trhdi' 



Tr 



F^\t,E)T(E)F^(t,E)-^ (T r (m;E) + Z°(m; E)) 



f(E) 



F^(t,E) ( U -^^f{E) + %Tm>{E)\ F^(t,E)^ (S7(m;£0 + Z a (m;E)) 



eV ar COS Ujt 



(All) 



ieh 



J 27rfi| 



Tr 



dF^(t,E) dF<-°M(t,E) 



dt 



-T(E) 



T m (E) 



dF^(t,E) dF<-°M(t,E) 



+Tr 



dE " iy ' dE 
dF(°\t,E) fdT(E) 



'T(E) 
df. 



dt 



T m (E) 



f(E) 



m , 0E f(E) + -^T(E)) F^(t,E)T m (E) (A12) 



0l r(E) ) OE^E)^ 
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h f dE^ fdF<®(t,E) (dY m ,{E) 



2irh' 



-Tr 



df- 



dE 

(F {0) (t,E)eV ac cosut 



oe m + ^w 



d_ 

dt 



(^F WA (t,E)eV ac costotj 



T m (E) 



-^E~ m + dE Tm ' {E) J dE Tm{E) (■ 



r 



(A13) 



The generalized injectiviti) 43 ^ 76 \ which is the density quence of the oscillating electrochemical potential in it, 
of particles on site i injected from the lead m' as a conse- is 



dn{i,Yuj,m r ) ie 



dV a , 



dE 



(f(E) - f(E + M) & A™, (E + hu- E)Ql - f(E)G 2 A r ^ (E + fiw; E)G 



(A14) 



2. Internal perturbation 



bias applied to the system. 



Now we will present the results for the response to an 
internal potential. Again, we assume that there is no DC 



Current response linear in V ac is given by 



dl m (lu) 

dVaa 



J dE { 



(f(E) - f(E + hu)) Tr \a%(E; E + hu)Q 2 WQ}\ - f(E)Tr [A£(£; E + tuv)Q 2 WQ ] 



+f(E + hu)Tr \a™(E; E + hiu)G\WGl 



(A15) 



The DC component of the current caused by the rec- as 
tification effect can be written (to feading order in V ac ) 



d 2 I m (0u>) 



= ^JdETv [(f(E) - f(E + fvu)) g Wg 2 T(E + hu)Q\W 'Q^ m (E) 

- (f(E -hu)- f(E)) g wg- 2 T(E - hw)gl 2 wglr m (E) 



(A16) 



Simifarly, the 2 nd harmonics of the current is given by 



d 2 I m (2to) _ e 3 



J dETr[(f(E) - f(E + hu)) G 2 W Q^_ 2 AZ{E -hu;E + hu) 



dV 2 c 2h 

(f(E -hu)- f(E)) G 2 WGoWgi 2 A%(E -hu;E + hu) + f(E + hu)G X 2 W G\W Gl 2 A a £{E -hu;E + hu) ( A17 ) 

-f(E - huj)G 2 WGoWG-2AZ{E - huj;E + hu) . 
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The response function to perturbing the onsite poten- 
tial on site i is the (generalized) emissiviti) 43 ^ 76 ^ and it 



reads (note, that it differs from the original definition by 
Biittiker, where one has to multiply it by l/(ieu)) 



dl m (lu) e 2 



h 



I 



dE 



(f(E) - f(E + M) GlA%(E; E + tuo)G 2 - f(E)g AZ(E; E + tuv)G 2 

+f(E + hio)glA%{E; E + hu3)Q\ . 



(A18) 



